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1.  INTRODUCTION 


Several  models  have  been  proposed  for  the  spall  process  that  accompanies  most 
nuclear  explosions.  The  models  are  very  difficult  to  validate  directly  because  (1)  the 
features  of  the  process  that  can  be  measured,  such  as  vertical  motion  at  the  surface,  are 
observed  at  only  a  few  locations  above  the  shot,  and  (2)  there  are  few  observations  of 
features  such  as  the  depth  and  area  of  the  spall  zone  and  the  distribution  of  motions 
within  the  zone.  Two-dimensional  nonlinear  simulations  of  the  explosion,  which 
include  the  physics  of  the  free-surface  interactions,  provide  an  opportunity  to  examine 
spall,  and  to  find  simple  models  to  represent  it.  We  present  a  linear  model  for  the 
spall  process  and  find  its  parameters  by  comparing  with  the  seismic  waves  from  explo¬ 
sion  simulations. 

Using  a  form  of  the  elastodynamic  representation  theorem,  we  have  calculated  the 
short-period  body  waves  emanating  from  the  source  zone  of  the  numerical  simulations. 
The  implications  for  teleseismic  magnitude  measurements  were  discussed  in  Day, 
et  al.  (1986)  and  McLaughlin,  et  al.  (1988).  By  comparing  these  P  and  SV 
wavefields  with  one-dimensional  nonlinear  calculations  (which  do  not  include  the  non¬ 
linear  effects  of  the  free-surface),  we  can  isolate  that  part  of  the  wave  field  due  to  the 
nonlinear  interaction  with  the  free-surface.  We  find  that  a  simple  tension  crack  with 
an  opening  that  propagates  with  the  pP  arrival  from  the  explosion  fits  P-waves  from 
simulations  of  both  Pahute  Mesa  and  Shagan  River  tests.  However,  the  model  which 
fits  the  P-wave  radiation  generates  SV-waves  that  are  too  small,  which  indicates  that 
there  are  additional  sources  of  shear  waves  in  the  finite  difference  simulations. 
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The  tension  crack  model,  proposed  in  its  original  form  by  Day,  et  al.  (1983),  is  a 
physical  model  whose  Darameters  can  be  compared  directly  to  field  observations.  The 
parameters  of  the  model  are  spall  depth  and  area,  detachment  velocity,  and  momentum. 
The  parameters  inferred  from  the  simulations  are  consistent  with  published  estimates  of 
these  values.  The  time  dependence  of  the  far-field  waveforms  is  a  natural  consequence 
of  the  model  and  requires  no  ad  hoc  choice  of  time  history  to  include  the  effects  of 
source  finiteness,  the  importance  of  which  has  been  pointed  out  by  Stump  (1985). 
Although  there  is  a  trade-off  between  the  parameters  of  the  model,  the  amplitudes  and 
waveforms  from  the  two-dimensional  (2D)  simulations  tightly  constrain  the  set  of 
parameters  which  fit  the  simulations. 

The  form  for  the  linear  source  representation  of  the  spall  model  is  such  that  it  can 
easily  be  added  to  a  one-dimensional  (ID)  explosion  source  to  compute  regional  and 
teleseismic  synthetic  seismograms.  In  a  companion  paper,  McLaughlin,  et  al.  (1990), 
use  the  model  to  compute  the  effects  of  spall  on  synthetic  regional  explosion  seismo¬ 
grams.  It  is  found  that  the  spall  contribution  to  the  Lg  signal  is  comparable  to  the 
direct  explosion  contribution. 

In  the  following,  we  describe  the  tension  crack  model  and  show  how  it  has  been 
implemented.  We  then  compare  the  body  waves  from  the  model  with  those  from  the 
numerical  simulations  and  compare  the  parameters  of  the  model  with  field  observa¬ 
tions.  Our  conclusions  are  presented  in  the  last  section. 
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2.  THEORY 


2.1.  Tension  Crack  in  a  Layered  Medium 

Our  model  for  the  spall  process  is  a  horizontal  circular  tension  crack,  shown 
schematically  in  Figure  1.  The  crack  lies  above  the  explosion  and  opens  when  the 
free-surface  rarefaction  from  the  explosion  encounters  the  crack.  The  material  above 
the  crack  is  assumed  to  fly  straight  up  and  return  under  the  influence  of  gravity  alone. 

We  have  extended  the  circular  tension  crack  model  of  Day,  et  al.  (1983)  to 
include  a  time  dependence  which  describes  the  detachment  and  slap-down  of  the  spall 
volume  as  well  as  the  finite  size  (radius)  of  the  crack.  The  spall  volume  in  the  model 
is  defined  as  the  cylinder  whose  radius  is  that  of  the  tension  crack  and  which  extends 
from  the  crack  to  the  surface.  Since  we  wish  to  compare  the  spall  model  with  2D 
simulations  which  were  done  in  a  layered  medium,  we  developed  the  formalism  for 
computing  the  far- field  body  waves  emanating  from  a  layered  medium.  The  parame¬ 
ters  of  the  model  are  the  crack  radius  and  depth,  and  the  distribution  of  detachment 
velocities  over  the  crack. 

In  this  section,  we  present  the  complete  body  wave  formalism.  We  then  make 
some  simplifying  assumptions  to  make  important  features  of  the  time  history  and  spec¬ 
trum  apparent. 

We  begin  with  the  representation  theorem  from  Aki  and  Richards  (1980),  which 
describes  the  displacement  field  due  to  a  discontinuity  on  a  surface  Z  : 

unii,t)  =  jj  mpq (£,/ )  *  Gnp^Lz-tt)  dZ 

where 
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Ground 

Surface 


X  is  the  observer’s  position, 

^  is  a  position  on  E, 

is  the  Green’s  function  giving  the  component  of  displacement  in 
the  xn  direction  at  x  due  to  a  force  in  the  t,p  direction, 

Gnp^ix-^, t)  is  the  gradient  of  G^ix-^t)  in  the  %q  direction, 
mpq  )  is  lhe  moment  tensor  density  for  the  discontinuity,  and 
*  denotes  temporal  convolution. 

For  an  isotropic  medium,  the  moment  tensor  density  is  given  by 


=  A.V*  [«*(£, OJ6w  +  p. 


vp[m<,(£,0!  +  Vq[up&t)] 


(2) 


where 


X  and  |i  are  the  Lame’  parameters  (at  £)  for  the  medium, 
lu  ]  is  the  slip  (crack  opening)  along  the  discontinuity, 
y  is  the  unit  normal  to  the  discontinuity,  and 
8pq  is  the  Kronecker  delta. 

(Repeated  indeces  indicate  summation  over  three  components).  For  a  horizontal  ten¬ 
sion  crack  with  opening  s  (^j,^2.t). 


¥7  =  [0,0,1] 

UfJr  =  [0,0,5(^2,O]  , 

and  the  moment  tensor  density  is 


(3) 


m 


X  0  0 

OX  0 

0  0  X+2p 


The  geometry  of  the  problem  is  shown  in  Figure  2. 
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f'i-iirc  2  Ciconutrv  lor  a  tension  crack  in  a  ’avered  medium. 
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We  now  specify  that  the  horizontal  crack  is  circular  with  radius  a  and  that  the 
discontinuity  varies  only  with  distance  from  its  center  (axially  symmetric).  Invoking 
cylindrical  symmetry  and  converting  to  cylindrical  coordinates,  the  displacement  field 
is 

2n  a 

Unu,t)  =//[>-  —  dra(roGnr0)  +  (X+2n)  aZoGWo]  *  s(r0,t)  r 0  dr0  d<p  .  (4) 

o  o  ro 

Here,  r0,<]),zo  are  the  source  cylindrical  coordinates 

4l  =  r o  cos(4») 
t,2  =  r0  sin(<J)> 

$3  =  z0 

Our  objective  is  to  compute  the  body  waves  leaving  a  layered  stack.  We  use  the 
formalism  of  Fuchs  (1966)  and  Bache  and  Harkrider  (1976).  The  Fourier  transforms 
of  the  P  wave  Green’s  functions  in  the  notation  of  Fuchs  are  given  by 

GtfCs.co)  =  Jf  Kfj  e-ikr^  kdk  (5) 

zn  0 

Here,  k  is  the  wavenumber  for  phase  velocity  c,  dJ  is  the  downgoing  wave  potential 
amplitude  due  to  a  force  in  the  j  direction,  F  is  the  vertical  distance  from  the  bottom 
of  the  stack  to  the  observer, 

ran  =  [(c/a„)2-l]1/2 

and  ka  =  o)/an  ,  where  a„  is  the  P  wave  speed  in  the  half-space  below  the  stack.  The 
details  on  the  computation  of  A;  can  be  found  in  Fuchs’  work.  The  kernels  Kfj  are 

Kzz  =  -ikr^JoikD  (6) 

K?z  =-iU](kF) 
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=-ikranJx{kr) 

K?r  =  ±k[J0m-J2(kr)] 

where  Jv  is  a  Bessel  function  of  the  first  kind  of  order  v.  The  coordinate 
F  =  I  r-rQ  I,  where  r  is  measured  from  the  center  of  the  crack  to  the  observer. 

The  SV  Green’s  functions  are 


CO)  =  -2(P„ Ic )2  f  o'  K&  e~ikr^T  kdk  ,  (7) 

o 

where  the  SV  potential  amplitudes  are 

K?zv  =  J0m  (8) 

K?2v  =  -irpnJx(kn 
K™ = jxm 

=  Yr^Jo(krhJ2m] 

Substituting  (5)  and  (6)  into  (4),  we  have  for  the  vertical  component  of  the  P 

wave 

2  it  a  » 

ufo, t)  =  — J  J  J  fl(*,F)  S(o,r0)  k2  r 0  cfr0  <«:  (9) 

ooo 

where 

B(k,r)  -  X  Ar°— [/tr^o^O+G-ro^Vi^n]  +  (*,+2*1)  dZoAz°J0(kr) 
ro 

Following  Fuchs  (1966),  we  evaluate  (9)  using  the  saddle  point  approximation,  which 
is  valid  for 

R  =  VzV2  »  yjz$+r$  ,  (10) 

That  is,  it  is  good  for  observer  distances  much  greater  than  the  source  radius  and 

depth.  Fuchs  shows  that  integrals  like  (9)  can  be  solved  using  the  relation 
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(11) 


J7(*,e» 


-ikja-1 


Jm(kr)dk  =/m+1  f  (k, 0))  e 


-ika-R 


ran  R  • 


Applying  this  to  Equation  (9),  we  have 


2n  a 


uPU,t)  =  -r-T  J  f  «  “  S(«Vo)  r0  dr0  d< j) 

Z7U(  0  jq 


(12) 


where 


/^(to,*)  =  X  Ar°k  +  (K+2[l)  3ZoAZo  , 

Cp  =  -cos 

and  QP  is  the  P  wave  take-off  angle.  The  wavenumber  k  is  evaluated  at  the  takeoff 
angle;  k=(ti/c ,  c=a„/sin(0/>).  To  evaluate  the  azimuthal  integral  in  (12),  we  invoke 
the  far-field  condition  (10),  and  approximate  the  slant  range  R  with  the  lower  order 
terms  of  its  Taylor  series. 


where 


Refering  to  Figure  2,  we  see  that 


R  =  r-~r  o 
r 


f  =  ylz^+r2 


—  =  sin0/>  cost})  , 
r 


so  that  (12)  becomes 


Cp 

vP, 


“z(x,0  =  I^rF  i<0,k)e 


2  7i  a 

ifcaaf  f  f  _  — iftrocos<t> 


■"  n e 
0  0 

Using  the  integral  representation  of  the  Bessel  function 


S((0,r0)  r0  dr0  dty  (13) 


wfll.O  =  ~  Fp(u),k)  e  lka*r  f  S (co,r 0)  J0(kr0)  r0  dr0  . 
K  o 

Applying  the  same  steps  to  the  SV  Green’s  functions  (7),  we  find 

csv  a 

u?V(£,t)  =  Fsv(w,k)  e~lk^r  f  S((0,rQ)  J0(krQ)  r0drQ 

K  rt 


where 


Fsv((a,k)  =  X  co r°k  +  (A.+2n)  aZocoz°  , 

C™  =  sin(20w)  . 

The  radial  integral  in  (14)  and  (15)  can  be  evaluted  numerically,  or  solved  analyt¬ 
ically  for  certain  algebraic  forms  of  S  (co,r0).  Some  examples  follow.  First  define 


P(co)  =  |  S((o,r0)J0{kr0)r0drQ  , 
o 


then  for  constant  crack  opening 


we  find 


5(co,r0)  =  S0(G>)  , 


P((o)  =  jJ1{ka)S0{(o)  . 

If  instead  the  dependence  of  crack  opening  on  radius  is  quadratic, 


5(co,r0)  =  50(co) 


ia2-r$) 


(co)  =  2k~2  Jiika)  50(co)  . 
Finally,  for  a  Gaussian  dependence, 


S  (co,r  0)  =  S0(o))e 


~(ro/a  )2 


/,(o»  =  ^  a2e'(fcj/2)2  50(co)  . 

We  note  that  the  effective  far-field  source  function  under  these  assumptions  is  the 


Hankel  transform  of  the  distribution  of  displacement  on  the  tension  crack.  We  shall 
refer  to  these  means  of  calculating  P(u>)  in  Section  3.  We  discuss  next  the 
specification  of  the  time  history  S(co,r0). 

2.2.  Time  History 

The  model  presented  here  is  based  on  a  tension  crack  which  opens  due  to  a  ten¬ 
sion  wave  from  the  free  surface.  We  assume  that  the  material  then  simply  behaves 
ballistically,  and  that  no  restoring  forces  in  the  medium  add  to  bring  the  material  back 
down  to  its  original  position.  The  displacement  at  each  point  on  the  tension  crack  then 
has  the  time  history 

_  •  V (r°’'o)T  "  \g *  f<*  0<  x  ^  2v (r0,l0)/g  (17) 

5  ^  ,r°^  “0  for  T  >  2v  (r0,tQ)/g 

b 

where  v(r0,r0)  is  the  detachment  velocity,  x=t~t0,  t0  is  the  time  of  initiation  at  r0  and 
g  is  the  acceleration  due  to  gravity.  The  form  of  v(r0,t0 )  can  be  regarded  as  a  param¬ 
eter  of  the  model,  and  we  have  experimented  with  various  forms  of  v(r0,t 0)  while 
fitting  the  model  to  numerical  simulations.  These  forms  included  (1)  analytic  distribu¬ 
tions  of  initial  velocities  which  initiate  at  the  same  instant  across  the  crack  ( tQ  is  con¬ 
stant),  (2)  random  distributions,  and  (3)  distributions  which  propagate  outward  with  the 
pP  phase  from  the  explosion.  The  first  form  allows  P  (to)  to  be  evaluated  analytically, 
as  shown  above  in  Section  2.1.  P( a)  must  be  evaluated  numerically  for  the  third 
form.  In  this  case,  we  set 

to(r0)  =  yl{2zs+zexp)2+r$  /  a 

and  v(r0)  to  some  algebraic  function  of  r0  .  As  we  shall  see  in  the  next  section,  the 
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analytic  first  form  was  adequate  for  the  Shagan  River  simulations,  while  the  third 
numerical  form  was  required  to  fit  the  Pahute  Mesa  simulations. 

Another  approach  is  to  view  the  spall  opening  as  a  stochastic  process  where  one 
has  only  probabilistic  estimates  of  the  process.  This  allows  us  to  make  general  state¬ 
ments  about  the  far-field  radiation.  We  assume  that  the  detachment  velocities  at  each 
point  on  the  crack  are  distributed  randomly  between  minimum  and  maximum  velocities 
v !  and  v2  according  to  a  probability  density  function  <j>(v).  We  make  the  bounds  expli¬ 


cit  by  writing 


where 


<()(v  )=/  (v )  [H  (v  -v  {)-H  (v  -v  2)]  . 


\  f(v)dv  =1. 


The  estimate  of  the  separation  history  is  then 


5(0=  J  t  (y-vg)  <(>(v)  dv 


where 


vs  =  TS< 


Introducing  the  moments  of  the  distribution 


Mj(a,b)  =  j  vj  f  (v)  dv  , 


the  separation  history  can  be  written  as 


5(0  =  1  fA/ ,  (v  i  ,v  2) — \AH  (0,/ , )  + 
t  [M  j  (vg  ,v  2)~Vg  M  0(vg  ,v  2)] AH  (t  j  ,t  2) 


where 
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and 


A H(tl,t2)  =  H(t~tl)-H(t-t2) 

is  a  boxcar’  turning  on  at  tj  and  off  at  t2.  We  have  included  the  details  of  the  open¬ 
ing  history  because  as  we  shall  see  in  the  next  section,  the  displacement  for  far-field 
waves  is  proportional  to  the  second  derivative  of  S(t),  which  is  given  by 

=  Mi(v„v2)  8(f)  -  gAH(0,tl)  +  (19) 

at 

l jvgf  (vg  ygM  0(vg  ,v2)]A H  (tht2) 

We  note  that  for  any  distribution  <b(v),  the  high  frequency  response  will  be  dominated 
by  the  8  function  term,  which  is  proportional  to  the  first  moment,  or  the  mean  velocity 
of  the  distribution.  This  delta  function  term  is  the  detachment  phase  which  is  followed 
by  a  free-fall  dwell,  which  is  the  second  term  in  the  equation.  The  last  term  is  slap 
down.  We  will  discuss  this  equation  and  other  features  of  the  model  in  more  detail  in 
the  following  sub-section. 

2.3.  Features  of  the  Model 

In  the  following,  we  make  some  simplifying  assumtions  for  the  purpose  of  illus¬ 
trating  the  features  of  the  model.  First,  we  assume  that  the  earth  model  is  a  uniform 
half-space.  In  this  case,  the  kernels  Fp  and  Fsv  become 

«■>.*)  -  [  <l-2p2p-2)  (^V^)  +  2P2ptla*,pe"“'  ]  (19) 

and 
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Fsv(a>,k)  =  [  2pr)asP  (e^-R^e^)  -  (\-2p2^2)  ]  (20) 

where  p-\lc  is  the  slowness,  and  ti„  =  a ~2-p2.  The  travel  times  for  the  P,  pP,  sP,  S, 
pS  and  sS  phases  are  denoted  by  tP ,  tpP ,  tsP ,  ts ,  tpS  and  ts$ ,  respectively.  The  vari¬ 
ables  RpP ,  RsP ,  RpS ,  and  RsS  are  the  corresponding  free  surface  reflection  coefficents 
(for  potentials),  given  by 

Rpp  =  [4TlaTlpp2  -  (p-2-2 P2)2]/  Dr  , 

Rsp  =-4np(r2-2p2)/^  , 

R„s  =4Tla(p-2-2p2)/D^  , 

with 

RsS  =  Rpp  • 

Dr  =  4tiaTipp2  +  ((T2-2p2)2 

Long  Period  Limits 

We  note  first  that  the  DC  values  of  Fp  and  Fsv  approach  zero  at  least  as  fast  as 
co2.  In  addition  to  the  factor  of  co  in  (19)  and  (20),  the  quantities 

d+2p2P"2)  (1  +RpP)  +  2$2pT]aRsP 
and 

2^a(l-RpS)-(l-2p2^2)RpS  . 

(the  limits  of  the  functions  in  brackets  in  (19)  and  (20)  as  co— >0),  are  identically  zero. 
As  pointed  out  by  Day  (1983),  the  spall  process  should  contribute  no  net  force  or 
momentum. 
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Small  Take-off  Angle  Approximation 

Next,  we  restrict  our  attention  to  steep  take-off  angles  (  small  8^  and  ®sv  )• 
Expanding  the  reflection  coefficients  to  second  order  in  take-off  angle,  we  find 

RpP  =  -1  +  (-j^-)30^  , 

Rsp  =  > 

Rg  =  -l  +  • 

Expanding  the  remaining  quantities  in  (19)  and  (20),  we  find  to  0(02y)  that 
' 

Fp  =  —  *  [l-2(-^-)26p]  (e~i<atp-e~iatpF)  -  8(-^)302(e“i<“''-e'ito,*p)  -  (21) 

a  a  a 

. 

Fsv  =  9S„  (e-‘““+e-i"*-2e-i“'s)  (22) 

Axial  symmetry  demands  that  the  SV  motion  be  zero  as  — >0,  as  seen  in  (22).  For 

small  angles,  the  P  wave  motion  is  dominated  by  the  P  and  pP  phases,  which  for  small 
depths  of  the  tension  crack,  appear  in  (21)  as  a  numerical  first  difference  operator. 

Approximation  for  Periods  Long  Compared  to  Travel  Time  to  the  Surface 

The  values  of  the  depth  of  the  tension  crack  that  we  infer  from  the  numerical 
simulations  in  Section  3  are  less  than  200  m.  The  P-pP  travel  time  for  the  Pahute 
Mesa  model  is  about  0.25  seconds,  and  less  than  one-half  that  for  the  Shagan  River 
model.  So  for  frequencies  much  less  than  about  4  Hz,  oi(tp-tpP)c\,  and 

e~iatr-e-iv*rr  =  2iw:r\azs  e~i(s*r  .  (23) 

Then  (21)  becomes  (to  0*A  order  in  QP  ), 
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(24) 


a 

Recall  Equation  (14),  for  the  P  wave  displacement. 


The  motion  along  the  ray  is 


approximated  by 


^(0))  =  -^Le"'^/>(co)  .  (25) 

y  crR 

Assuming  that  crack  opening  is  constant,  and  taking  the  inverse  Fourier  transform  of 
(25),  we  obtain 


u£v(0  =  ~T~  [(ap)2~t2]m  *  SQ(t-tP)  .  (26) 

pzvSR 

Finally,  we  assume  that  the  periods  of  interest  are  also  longer  than  ap ,  the  time  it 
takes  for  a  wave  to  travel  across  the  tension  crack.  Typical  teleseismic  values  of  p  are 
less  than  0.1  sec/km  and  values  of  a  are  around  1  km,  so  pa  <0.1  sec.  In  this  case,  we 
can  write 

P~2  [(ap)2-t2]m  -  y7ta2  5(?)  , 

(26)  becomes 


uray(0  =  2(\+2[L)R  ^  'p)  ’ 


(27) 


where  mspau  =pnzs  a 2  is  the  mass  of  the  material  lifted  by  the  spall  process.  The  P- 
wave  amplitude  is  then  proportional  to  the  spall  mass,  and  the  opening  acceleration  S . 
Recall  that  in  Section  2.2,  we  proscribed  the  opening  history  for  a  distribution  of 
detachment  velocities.  We  assume  the  distribution  is  uniform  (f  (v)=l.)  between  velo¬ 
cities  vt  and  v2,  evaluate  (19),  and  substitute  the  result  in  (27): 


Uray(t)  = 


m 


spall 


2(k+2\x)R 


[v  5(r’)  -  gAH(0,t[)  + 


(28) 
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— S —  (.-V  2+-J  gt' )  AH  (t[,t  2)] 
v2-vt  4 

where  t=t-tP  and  V=(vx+v2)/2.  Equation  (28)  is  shown  schematically  in  Figure  3. 
The  inital  8  function  term  is  proportional  to  mspallv  ,  the  mean  spall  momentum.  The 
dwell  term  is  mspaI[g  ,  the  gravitational  restoring  force.  The  duration  of  the  dwell  is 
fj,  and  that  of  the  slap  down  is  tx-t2. 

The  far-field  SV  displacement  is  approximated  by 

U%(t)  =  ^  esv  SCO  *  m-ts)+8(t-tsS)-28(t-tpS)}  (29) 

Spectral  Characteristics 

The  spectrum  of  the  body  waves  from  the  tension  crack  model  is  a  narrowband 
signal.  To  see  this,  consider  Equation  (25).  As  co-*0,  P  (to)  approaches  a  constant  pro¬ 
portional  to  the  area  of  the  crack,  so  the  far-field  displacement  approaches  zero  as  to2. 
This  is  the  case  for  any  type  of  dependence  of  opening  on  r0.  The  frequency  limit 
depends  on  the  details  of  5(co,r0),  but  we  can  make  some  general  statements.  Using 
the  time  dependence  discussed  above  in  Equation  (19),  it  can  be  shown  that  at  high 
frequencies,  S  (to)=M  ,(v  j.Vj),  the  mean  detachment  velocity.  Then 

a 

P(co)  =  Af  J  J  S  (r  0)  r  0  dr0  . 

0 

For  constant  opening  (Section  2.1),  upon  using  the  asymptotic  approximation  for  the 
Bessel  function, 

MraV(co)  =  (o-3/2  -—-e~,w,pP(03)  cos(ka- 3n/4)  ,  (25) 

*  a  2R 

which  rolls  off  at  a  rate  of  to_3/2  at  high  frequencies.  However,  for  a  crack  opening 
distribution  which  goes  to  zero  at  the  edge  of  the  crack,  the  roll  off  is  faster.  Consider 
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/ 


Figure  3.  Schematic  drawing  of  the  detachment  acceleration  history  which  is 
proportional  to  the  far-field  P-wave  displacement. 
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a  distribution 


S0(r0)  =  1  -  • 

a 

where  n  is  a  positive  integer.  With  a  little  manipulation, 

ka 

-i-  , 

(*a)n  £  o 

where  £=&r0.  It  can  be  shown  that  at  high  frequencies,  this  expression  is  proprotional 
to  k~2,  so  that  the  displacement  rolls  off  at  co-5/2.  Thus,  any  radial  dependence  which 
goes  to  zero  as  fast  as  a  power  of  r0,  will  have  a  roll  off  of  gT5/2.  The  spectrum  is 
shown  schematically  in  Figure  4.  Since  explosion  models  typically  decay  as  to~2  at 
high  frequencies,  and  are  flat  at  low  frequencies,  the  tension  crack  will  be  a  more  nar¬ 
rowband  signal  than  the  explosion  by  itself.  Thus,  if  the  spall  contribution  is  compar¬ 
able  to  that  of  the  explosion,  it  will  be  so  only  in  a  narrowband. 
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Log  Frequency 


I  imue  4.  Schematic  drawing  of  far-held  P-vvavc  displacement  spectrum. 
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3.  COMPARISON  WITH  NUMERICAL  SIMULATIONS 


The  objective  of  this  study  was  to  derive  a  simple,  linear  physical  model  which 
represents  the  spall  process.  Our  strategy  has  been  to  isolate  source  interaction  effects 
by  comparing  the  wavefields  from  2D  calculations,  «2D>  those  from  ID  calcula¬ 
tions  done,  u  w ,  in  the  same  source  materials  (but  in  an  inifinite  medium).  Our  work¬ 
ing  assumption  is  that  the  difference  in  the  wavefields  represents  the  free  surface 
interaction  effects,  as  sketched  in  Figure  5.  In  fact,  the  procedure  has  been  to  find  the 
difference  in  the  displacements  uw-u 1D ,  and  vary  the  parameters  of  the  model 
(described  in  the  previous  section)  to  fit  the  displacements  uspall  .  Alternatively,  we 
could  have  found  uspall  which  satisfies  uw+uspaU  =  «2D-  We  found  the  first  approach 
to  be  more  direct. 

To  find  the  ground  motions  from  the  2D  calculations,  we  used  the  methods 
described  in  Rodi,  et  al.  (1978),  Bache,  et  al.  (1982),  and  Day,  et  ai.  ,  (1983, 
1986).  The  2D  simulations  themselves  are  described  in  Day,  et  al.  (1986).  The  elas¬ 
tic  properties  for  the  Pahute  Mesa  and  Shagan  River  simulations  are  shown  in  Table 
3.1.  The  method  for  computing  the  body  waves  entails  evaluating  a  form  of  the  elas- 
todynamic  representation  theorem,  which  gives  the  motions  in  terms  of  spatial  and 
temporal  convolutions  of  displacements  and  stresses  monitored  on  a  surface  surround¬ 
ing  the  nonlinear  zone  with  Green’s  functions  and  their  spatial  gradients. 
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Spall  Zone 
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Figure  5.  Sketch  of  the  P-waves  from  the  2D  nonlinear  simulations  (left)  and  ID  cases  (right)  in  which  an  RDP  (from  a 
ID  nonlinear  simulation)  is  embedded  in  an  elastic  medium. 


Table  3.1 

Elastic  Structure  for  die  Shagan  River  Simulations 

T°yer 

Compressional 

Shear 

Thickness 

Velocity 

Velocity 

HI 

(m) 

(m/sec) 

(m/sec) 

oo 

5018 

2789 

2700 

Elastic  Structure  for  the  Pahute  Mesa  Simulations 

Layer 

Compressional 

Shear 

Thickness 

Velocity 

Velocity 

HI 

(m) 

(m/sec) 

(m/sec) 

112.5 

1208 

661.4 

1600 

457.5 

2025 

1109 

OO 

2887 

1581 

3.1.  Shagan  River  Simulations 

We  begin  with  the  Shagan  River  simulations,  which  were  all  done  at  a  yield  of 
125  KT  and  at  three  depths  of  burial  (DOB):  an  over  buried  depth  (980  m),  an 
optimally  buried  depth  (680  m),  a  depth  just  below  cratering  (no  ejecta)  (300  m),  and  a 
depth  which  causes  cratering  (200  m).  The  results  for  P  waves  were  presented  in  Day, 
et  al.  (1986),  and  we  include  them  here  for  purposes  of  comparing  with  the  spall 
model.  As  can  be  seen  in  Table  3.1,  the  earth  model  is  a  half-space,  with  material 
properties  representing  typical  values  for  the  high  velocity  near  surface  media  at  the 
site.  The  far-field  P  and  SV  displacements  are  shown  for  three  take-off  angles  in  Fig¬ 
ures  6  to  11.  The  take-off  angles,  10°,  20°  and  30°,  are  representative  of  teleseismic, 
intermediate  and  regional  slownesses  (phase  velocities)  for  this  source  structure.  The 
motions  from  the  2D  simulations  are  overlain  with  ID  simulations  in  which  the  non¬ 
linear  properties  at  the  source  are  the  same  as  the  2D  values.  The  ID  source  (RDP) 
was  inferred,  and  the  body  waves  were  calculated  for  a  linear  elastic  medium.  Thus, 
for  the  Shagan  River  half-space  model,  the  ID  waveforms  have  the  direct  P  phase,  as 


23 


Id,  dob-200  _  Id,  dob-300 
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Comparison  of  the  far-field  P-waves  from  the  ID  (solid  lines)  and  2D  (dashed  lines) 
Shagan  River  simulations  for  four  depths  of  burial  at  a  take-off  angle  of  30*. 
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Id,  dob-200  _  Id,  dob-300 


Comparison  of  the  far-field  S  V-waves  from  the  ID  (solid  lines)  and  2D  (daslrrd  lines) 
Shagan  River  simulations  for  four  depths  of  burial  at,  a  take-off  angle  of  20°. 


well  as  the  pP  and  pS  elastic  reflections.  These  signals  do  not  include  the  effects  of 
mantle  or  crustal  propagation,  anelastic  attenuation,  or  recording  instrumentation. 

In  each  case,  the  rise  time  and  amplitude  of  the  first  peak  in  the  P  waves  (Figures 
6-8)  is  very  nearly  the  same  for  the  ID  and  2D  calculations.  For  the  shallow  calcula¬ 
tions,  200  and  300  m,  the  ID  and  2D  waveforms  diverge  after  the  peak  and  are  quite 
different  at  later  times,  indicating  that  the  cratering  and  near-cratering  processes  are 
indeed  different  from  elastic  pP  reflections.  For  the  deeper  DOB’s,  680  and  980  m, 
the  apparent  pP  phase  on  the  2D  records  appears  to  be  diminished  and  delayed  relative 
to  the  ID  case. 

The  SV  waves  (Figures  9-11)  show  much  greater  differences  between  the  ID  and 
2D  cases.  The  2D  solutions  have  a  direct  S  wave,  due  to  vertical  asymmetries  in  the 
source,  which  is  not  in  the  ID  solutions.  The  main  peaks  in  the  2D  SV  waves  are 
larger  than  the  ID  peaks  at  all  four  DOB’s,  and  are  about  twice  as  big  for  300,  680 
and  980  depths.  The  duration  of  the  main  peaks  is  also  greater  for  the  2D  cases.  A 
large  negative  swing  occurs  on  the  2D  waveforms,  which  is  much  smaller  on  the  ID 
signals. 

As  discussed  above,  we  computed  the  2D- ID  waveforms  and  matched  them  to 
those  from  the  spall  model.  The  parameters  of  the  spall  model  which  most  closely 
matched  the  numerical  simulations  are: 
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Table  3.2 

Spall  Model  Parameters  for  the  Shagan  River  Simulations 

Depth 

1 

Crack 

Minimum 

Maximum 

of  Burial 

■ 

Radius 

Detachment 

Velocity 

Detachment 

Velocity 

zexp 

a 

vi 

v2 

(m) 

■m 

(m) 

( m/sec ) 

(m/sec) 

200 

100 

600 

1.1 

20.0 

300 

150 

600 

1.1 

15.0 

680 

200 

600 

1.1 

4.0 

980 

600 

1.1 

1.5 

The  theoretical  waveforms  were  found  to  be  weakly  dependent  on  the  form  of  the 
source  function  S(ay0)  in  this  case,  so  it  was  assumed  constant.  In  Figures  12 
through  15,  we  compare  the  P  wave  spall  model  waveforms  with  the  2D- ID 
waveforms  for  each  of  the  DOB’s.  In  general,  the  comparisons  are  quite  good,  both  in 
amplitude  and  shape.  The  parameter  that  varies  the  most  in  Table  3.2  is  the  maximum 
detachment  velocity  v2  .  Recall  that  in  Equation  (28),  Section  2.3,  v2  enters  the  solu¬ 
tion  in  two  ways.  First,  the  5  function  term,  which  causes  the  first  peak  in  the 
waveforms,  is  proportional  to  the  mean  detachment  velocity  v".  Second,  the  duration 
of  the  signal  (end  of  slap-down)  is  time  r2  =  2 v2/g.  As  the  DOB  decreases,  the  2D- 
1D  difference  waveforms  increase  in  amplitude  and  duration.  Thus,  as  the  DOB 
decreases,  v2  is  required  to  increase.  We  note  that,  even  in  the  case  where  cratering 
occured  (200  m),  the  tension  crack  model  provided  a  good  representation  of  the  P  radi¬ 
ation.  The  shapes  of  the  signals  for  the  980  m  case  agree  well.  For  the  680  m  case, 
the  agreement  is  good  except  for  an  additional  inflection  in  the  later  parts  of  the  2D 
waveforms  that  is  not  modeled  by  the  tension  crack. 
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The  peak  detachment  velocity  is  inferred  from  field  observations  from  ground 
motion  records  at  or  near  ground  zero.  We  can  do  an  analogous  measurement  by  exa¬ 
mining  the  vertical  velocity  in  the  finite  difference  simulation  at  the  free-surface  of  the 
grid.  We  find  that  the  values  of  v2  in  Table  2.2  agree  closely  with  the  ground  zero 
velocities  in  the  simulations. 

Using  the  parameters  in  Table  3.2,  we  computed  the  SV  waves  from  the  tension 
crack  model  and  compared  them  with  the  corresponding  2D-1D  SV  waves.  As  can  be 
seen  in  Figures  16  to  19,  the  comparisons  are  not  as  good  as  the  P-waves,  especially  at 
the  deeper  DOB’s,  where  the  amplitudes  predicted  by  the  model  are  too  small  by  fac¬ 
tors  of  two  to  three.  Attempts  to  find  a  set  of  parameters  which  would  improve  the 
SV  comparisons  without  degrading  the  P  comparisons  were  unsuccesful.  We  note  that 
in  all  cases,  the  tension  crack  models  derived  from  the  P-waves  underestimate  the  SV 
radiation. 

The  results  of  the  SV  comparisons  suggest  that  the  tension  crack  model  is  not 
rich  enough  in  SV  radiation.  We  therefore  tried  some  additions  to  the  model  to 
improve  this  situation.  The  additions,  a  compensated  linear  vector  dipole  and  a  distri¬ 
buted  double-couple,  are  pictured  in  Figure  20.  We  hypothesized  that  processes  tran¬ 
spired  in  the  simulations  that  could  be  represented  as  modifications  to  the  isotropic  part 
of  the  moment  tensor.  We  first  chose  to  represent  this  as  a  compensated  linear  vector 
dipole  or  CLVD  (  e.g.,  Knopoff  and  Randall,  1970).  The  CLVD  represents  a  process 
akin  to  squeezing  a  vertical  tube  of  toothpaste.  The  moment  density  for  the  CLVD  is 
given  by 
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Figure  20.  Sketch  of  alternative  sources  of  SV:  the  compensated  linear  vector  dipole  (left)  and  a  distributed  double¬ 
couple  (right). 
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This  moment  tensor  has  zero  trace.  Following  the  analysis  in  Section  2,  the  formulas 
for  the  displacements  (P  along  the  ray,  SV  normal  to  the  ray)  are  similar  to  Equations 
(14)  and  (15), 
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We  then  formed  the  sum 


^  ^  ^ crack  H CLVD 

where  McracJt  is  the  tension  crack  solution,  and  0<v<l.  We  assumed  the  time  and  spa¬ 
tial  dependence  of  S  (to ,r0)  was  the  same  for  both  the  tension  crack  and  CLVD.  Fig¬ 
ures  21  and  22  show  the  comparisons  for  the  P  and  SV  waves  for  the  DOB=980  case 
with  v  =  0  (all  CLVD,  no  tension  crack),  which  we  found  works  best  for  this  DOB. 
We  see  that  the  P  wave  solutions  are  changed  little  except  at  the  shallowest  take-off 
angle  (  30°  ),  where  the  second  peak  in  the  CLVD  solutions  is  now  too  large.  On  the 
other  hand,  the  comparisons  for  the  SV  waves  are  greatly  improved.  The  shape  and 
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Comparison  of  the  far-field  SV-wavcs  from  the  2D-1D  (solid  lines)  Shagan  River 
simulations  with  the  CLVD  model  (dashed  lines)  (with  no  tension  rrark  contribution)  at 
three  take-off  angles  for  a  depth  of  burial  of  980  m 
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Figure  22.  Comparison  of  the  far-field  P-waves  from  the  2D-1D  (solid  lines)  Shagaii  Ri 
simulations  with  the  CLVD  model  (dashed  lines)  (with  no  tension  rraek  contribution) 
three  take-off  angles  for  a  depth  of  lmrial  of  980  m 


amplitudes  agree  fairly  well.  The  initial  negative  part  of  the  2D- ID  signals  is  a  direct 
S  wave  from  the  explosive  source.  It  is  not  part  of  the  free  surface  interaction  which 
we  are  modeling  here.  For  the  DOB=680  case,  we  found  that  a  value  of  v  =  0.5 
(equal  parts  tension  crack  and  CLVD)  works  better  than  other  choices,  but  the  results 
are  not  as  good  as  for  DOB=980,  as  can  be  seen  in  Figures  23  and  24.  Although  the 
amplitudes  of  the  main  peak  are  matched  well,  the  apparent  period  and  change  with 
take-off  angle  of  the  following  trough  are  matched  poorly.  We  find  that  adding  the 
CLVD  source  to  the  tension  crack  for  the  shallow  DOB’s  (200  and  300)  worsens  the 
comparisons. 

Since  this  approach  showed  some  promise,  we  also  hypothesized  that  additional 
sources  of  shear  waves  may  be  emanating  from  vertical  cracks  which  are  the  boun¬ 
daries  of  a  cylinder  which  moves  above  the  horizontal  tension  crack.  That  is,  we  ima¬ 
gine  that  the  material  above  the  tension  crack  moves  up  and  down  as  a  unit,  shaped  as 
a  cylinder  whose  base  is  the  tension  crack.  The  relative  motion  along  the  vertical 
sides  would  act  as  a  distributed  shear  dislocation.  Without  showing  the  details,  we 
summarize  by  saying  that  after  including  this  source  in  the  computer  code,  its  contri¬ 
bution  to  the  DOB=680  and  980  SV  solutions  was  too  small  to  match  the  2D- ID 
differences. 

The  tension  crack  model  presenteu  here  matches  the  P  waves  due  to  free  surface 
interactions  in  these  simulations  very  well,  but  the  model  includes  only  part  of  the  SV 
wave  generation  mechanisms. 
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Comparison  of  2D- ID  far-field  SV-wavos  (solid  linos)  and  results  using  1/2  CLVD  and 
l/2  tension  crack  contributions  (dashed  lines)  at  three  take-off  angles  for  a  depth  of 
burial  of  080  m  from  the  Shagari  River  simulations. 


3.2.  Pahute  Mesa  Simulations 


The  2D  and  ID  P-wave  signals  are  overlain  in  Figures  25  through  27.  The  calcu¬ 
lations  were  done  at  DOB’s  of  200,  680  and  980  m  for  the  Pahute  Mesa  model,  but 
not  at  300  m.  The  take-off  angles  (3. 3°,  11.4°,  16.7°)  were  chosen  so  that  the  phase 
velocities  are  the  same  as  in  the  Shagan  River  simulations  above.  The  first  peaks  of 
the  ID  and  2D  signals  align  closely.  For  DOB’s  680  and  980,  the  apparent  pP  from 
the  2D  waveforms  is  lagged  about  0.2  seconds  relative  to  the  elastic  pP  in  the  ID 
waveforms.  The  2D  pP  amplitudes  are  attenuated  relative  to  ID  amplitudes  for  the 
shallower  take-off  angles  (11.4°  and  16.7°). 

As  with  the  Shagan  River  simulations,  we  fit  the  tension  crack  model  to  the  2D- 
1D  difference  time  series.  The  results  are  shown  in  the  comparison  plots  in  Figures  28 
through  30.  The  fits  of  the  model  to  the  2D- ID  signals  are  generally  good  for  the 
three  DOB’s  and  take-off  angles.  The  parameters  of  the  model  are 


Table  3.3 


Spall  Model  Parameters  for  the  Pahute  Mesa  Simulations 


Crack 

Crack 

Minimum 

Maximum 

IJBSllIfi 

Depth 

Radius 

Detachment 

Detachment 

■ 

Velocity 

Velocity 

'WBmS. 

ZS 

a 

vi 

v2 

(m) 

(m) 

(m/sec) 

(m/sec) 

150 

400 

1.1 

22 

200 

2200 

1.1 

2 

150 

2200 

1.1 

2 

There  are  several  differences  between  the  model  parameters  derived  here  and  those  for 
the  Shagan  River  simulations.  First,  in  the  Shagan  case,  we  found  that  the  solutions 
were  insensitive  to  the  choice  of  the  form  of  the  source  function  S  ((0,r0).  This  is  not 
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Comparison  of  thr  far-firbl  P-wavos  from  tlir  ID  (> 
Pahuto  Mrsa  simulations  for  tlirrr  depths  of  burial  at 


e  27.  Comparison  of  the  far-field  P-waves  from  the  ID  (solid  lines)  and  2D  (d 
Pahute  Mesa  simulations  for  three  depths  of  burial  at.  a  take-off  angle  of  16.7 
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Figure  29.  Comparison  of  the  far-Field  P-waves  from  the  2D-1D  (solid  lines)  Pahutc  Mesa  simulations 
with  the  tension  crack  model  (dashed  lines)  at  three  take-off  angles  for  a  depth  of  burial 
of  680  m. 


the  case  for  the  Pahute  Mesa  runs,  where  we  found  that  S((o,r0)  must  be  distributed 
over  the  crack  with  the  largest  values  near  the  center.  This  is  due  to  the  lower  wave 
speeds,  and  hence  phase  velocities,  in  the  Pahute  Mesa  model  (Table  3.1).  In  Figures 
28  through  30,  a  Gaussian  spatial  dependance  with  onset  lagged  by  the  pP  time  from 
the  explosion  (see  Section  2.1  and  2.2)  was  used.  The  solutions  were  not  sensitive  to 
the  details  of  the  spatial  dependance,  as  long  as  the  maximum  was  at  the  center  and 
went  to  zero  at  the  edge.  This  spatial  dependance  was  required  to  match  the  change  in 
2D- ID  signals  with  take-off  angle.  It  was  not  needed  for  the  Shagan  River  runs 
because  the  wave  speeds  were  about  twice  those  for  Pahute  Mesa.  For  example,  we 
show  the  tension  crack  signals  for  a  constant  S(a>,r0)  for  DOB=980  with  the  2D- ID 
signals  in  Figure  31.  It  can  be  seen  that  the  results  are  not  as  good  as  in  Figure  29 
where  the  Gaussian  spatial  dependance  was  used. 

Another  feature  of  the  Pahute  Mesa  simulations  that  differs  from  the  Shagan 
River  runs  is  the  large  change  in  crack  radius  between  the  case  for  DOB=200  and  the 
other  two  DOB’s.  This  is  presumably  due  to  the  shallow  layering  in  the  Pahute  Mesa 
earth  model.  Although  the  tension  crack  had  a  maximum  detachment  velocity  (22 
m/sec)  similar  to  Shagan,  it  occured  over  a  much  smaller  radius  (400  m). 

A  further  difference  between  the  models  for  the  two  sites  is  that  the  maximum 
detachment  velocities  for  Pahute  Mesa  were  the  same  for  DOB’s  680  and  980, 
whereas  the  velocity  decreased  from  DOB=680  to  980  in  the  Shagan  model.  This  is 
again  attributed  to  differences  in  the  shallow  nonlinear  properties. 
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We  computed  the  SV  waves  for  the  tension  crack  using  the  parameters  in  Table 
3.3,  and  compared  them  to  the  2D- ID  difference  waveforms.  They  are  shown  in  Fig¬ 
ure  32,  for  DOB=680  m.  As  with  the  Shagan  River  simulations,  the  2D- ID  signals 
are  larger  than  those  from  the  tension  crack.  Again  this  indicates  that  processes  are 
generating  SV  waves  in  the  simulations  that  are  not  included  in  the  tension  crack 
model. 


3.3.  Comparisons  with  Observations 

From  observations  of  surface  ground  motions  and  of  physical  manifestations  of 
Pahute  Mesa  tests,  several  authors  have  estimated  the  parameters  of  the  spall  process. 
In  Table  3.4,  we  compare  their  results  with  those  for  the  tension  crack  model  at  the 
optimal  DOB  (680). 


Table  3.4.  Spall  parameters  from  field  observations  and 
from  the  tension  crack  model  for  Pahute  Mesa 

Parameter 

Tension 

Crack 

Patton 

(1990) 

Sobel 

(1978) 

Viecelli 

(1973) 

■Bn 

Rawson 

(1988) 

Maximum 

Velocity 

(m/sec) 

2 

D 

8 

Radius 

(m) 

2200 

1750 

500 

500- 

2665 

Momentum 
(  xlO12  Nt-s) 

m 

9.2 

3.5 

0.58 

Mass 

(xlO12  Kg) 

5.3 

3.5 

1.2 

0.2 

3.0 

Depth 

(m) 

200 

110 

100- 

400 

The  authors  cited  in  the  table  typically  expressed  their  results  as  scaled  values,  in 
which  case  we  used  125  KT,  the  yield  in  the  simulations.  In  addition,  the  values  in 
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Table  3.4  are  the  average  values  given  by  the  authors.  The  mean  detachment  velocity 
(1.55  m/sec)  was  used  to  compute  the  tension  crack  momentum,  rather  than  the  max¬ 
imum  velocity.  The  momentum,  spall  radius  and  depth  from  the  tension  crack  model 
lie  within  the  spread  of  those  inferred  from  observations.  The  mass  is  at  the  high  end 
while  the  maximum  velocity  is  at  the  low  end. 

It  is  very  difficult  to  estimate  the  parameters  in  Table  3.4  because  of  a  paucity  of 
direct  observations.  There  are  typically  few  surface  ground  mo^on  sensors,  and  no 
sub-surface  observations  of  ground  motion  or  other  phenomena  such  as  cracking. 
Inference  of  the  spall  process  beyond  estimating  apparent  pP  amplitude  and  travel  time 
from  far-field  recordings  has  not  been  successful.  With  these  limitations  in  mind,  it 
can  be  sa»d  that  the  parameters  of  this  study  are  consistent  with  published  field  obser¬ 
vations. 
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4.  SUMMARY  AND  CONCLUSIONS 


We  show  that  a  simple  model  can  be  constructed  which  generates  the  far-field  P 
waves  of  a  two-dimensional  nonlinear  calculation  which  includes  the  effects  of  the 
free-surface  for  a  range  of  depths-of-burial.  The  model  which  fits  the  P-waves  from  the 
Shagan  River  simulations  underestimates  the  SV-waves.  The  model,  proposed  in  its 
original  form  by  Day,  et  al.  (1983),  is  based  on  using  a  tension  crack  which  opens 
due  to  the  tension  wave  from  the  free  surface.  The  material  over  the  crack  travels  up 
with  the  impulse  of  the  tension  wave  and  returns  nder  the  influence  of  gravity. 
Stump  (1985)  modified  the  original  model  to  include  the  effects  of  source  by  adding 
an  empirical  time  function  based  on  chemical  explosions.  The  modified  model  was 
used  by  Taylor  and  Randall  (1989)  to  model  regional  seismograms.  In  our  formalism 
presented  in  this  report,  we  include  the  effects  of  source  finiteness  and  those  of  crustal 
reverberations.  A  time  dependence  is  a  natural  consequence. 

The  parameters  of  the  model  compare  favorably  with  observations  based  on  field 
data.  In  addition,  McLaughlin,  et  al.  (1990),  use  the  model  to  compute  regional 
seismograms  and  find  that  in  the  Lg  bandwidth,  the  spall  contribution  should  be  com¬ 
parable  to  or  greater  thar  that  of  the  explosion  signal  alone.  Since  the  spall  model 
which  Ins  the  P-waves  in  our  study  generates  SV-waves  that  are  too  small,  the  spall 
model  may  lead  to  an  underestimate  of  the  Lg  from  spall.  Spall,  therefore,  appears  to 
be  a  very  significant  source  of  Lg. 
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